arXiv: 1504.0135 lv3 [physics.plasm-ph] 18 Aug 2015 


Gyrokinetic and kinetic particle-in-cell simulations of guide-field 
reconnection. I. Macroscopic effects of the electron flows 

P. A. Munoz, 1 ' a) D. Told , 2 ' 3 P. Kilian , 1 J. Buchner , 1 and F. Jenko 2 ’ 3 

14 Max-Planck-Institut fur Sonnensystemforschung, D-37077 Gottingen, Germany 
^ Max-Planck-Institut fur Plasmaphysik, D-857J8 Garching, Germany 
Department of Physics and Astronomy, University of California, Los Angeles, California 90095, 

USA 

(Dated: August 19, 2015) 

In this work, we compare gyrokinetic (GK) and fully kinetic Particle-in-Cell (PIC) simulations of magnetic 
reconnection in the limit of strong guide field. In particular, we analyze the limits of applicability of the 
GK plasma model compared to a fully kinetic description of force free current sheets for finite guide fields 
(b g ). Here we report the first part of an extended comparison, focusing on the macroscopic effects of the 
electron flows. For a low beta plasma (/3j = 0.01), it is shown that both plasma models develop magnetic 
reconnection with similar features in the secondary magnetic islands if a sufficiently high guide field ( b g > 30) 
is imposed in the kinetic PIC simulations. Outside of these regions, in the separatrices close to the X points, 
the convergence between both plasma descriptions is less restrictive ( b g > 5). Kinetic PIC simulations using 
guide fields b g < 30 reveal secondary magnetic islands with a core magnetic field and less energetic flows 
inside of them in comparison to the GK or kinetic PIC runs with stronger guide fields. We find that these 
processes are mostly due to an initial shear flow absent in the GK initialization and negligible in the kinetic 
PIC high guide field regime, in addition to fast outflows on the order of the ion thermal speed that violate 
the GK ordering. Since secondary magnetic islands appear after the reconnection peak time, a kinetic 
PIC/GK comparison is more accurate in the linear phase of magnetic reconnection. For a high beta plasma 
(/3j = 1.0) where reconnection rates and fluctuations levels are reduced, similar processes happen in the 
secondary magnetic islands in the fully kinetic description, but requiring much lower guide fields ( b g < 3). 
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I. INTRODUCTION 

Magnetic reconnection is a fundamental process in 
plasma physics that converts magnetic field energy ef¬ 
ficiently into plasma heating, kinetic energy of bulk flows 
and accelerated particles. It plays a role in a wide range 
of different environments, from fusion devices, planetary 
magnetospheres, and stellar coronae to extragalactic ac¬ 
cretion disks. 1-3 

Originally, magnetic reconnection was modeled with 
antiparallel magnetic fields in a plane sustained by a cur¬ 
rent sheet (CS). But that configuration is not common 
in nature, since the presence of an out-of-plane magnetic 
field is ubiquitous in both space and laboratory plasmas, 
such as in the solar corona and fusion devices, respec¬ 
tively. However, the properties of reconnection in the 
presence of a guide field (b g ) are still nowadays much 
less understood. 4 Thus, the study of stability of these 
CS with shear magnetic fields is of paramount impor¬ 
tance to understand the conversion of magnetic energy 
in these environments. 

In order to capture the kinetic effects in magnetic 
reconnection, fully kinetic particle-in-cell (PIC) codes 
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have become one of the preferred tools since some time 
ago, with the increasing availability of computational 
resources. 5 These codes allow the simulation of magnetic 
reconnection in fully collisionless Vlasov plasmas. How¬ 
ever, these fully kinetic PIC simulations (for brevity, from 
now on these ones will be simply referred to as “PIC sim¬ 
ulations”) can be computationally very demanding due 
to stability reasons and the requirement of keeping nu¬ 
merical noise as low as possible. One approach to solve 
this problem is using gyrokinetic (GK) theory, an ap¬ 
proximation to a Vlasov plasma (see Ref. 6 or Refs. 7-9 
for more recent reviews). In this plasma model, the idea 
is to decouple the fast gyrophase dependence from the 
dynamics, considering the motion of charged rings, thus 
reducing the phase space from six to five dimensions and 
removing dynamics on very small space-time scales. The 
gyrokinetic approach is suitable in strongly magnetized 
plasmas (equivalent to the limit of very large magnetic 
guide fields in magnetic reconnection). More precisely, 
the ion gyroradius pi has to be much smaller than the typ¬ 
ical length scale Lq of the background: e = pi/Lg <C 1, 
where e is the GK ordering parameter. In contrast, it is 
important to note that the typical length scale of vari¬ 
ation of the perpendicular perturbations, kj 1 , can be 
comparable to the ion gyroradius: k±pi ~ 1. The order¬ 
ing e 1 implies the restriction of the GK approach to 
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phenomena with frequencies lower than the ion cyclotron 
one: u>/Cl C i ~ 0(e), where w is a typical frequency in the 
system. The fluctuations in the distribution function and 
electromagnetic fields, with respect to their equilibrium 
values, are also of order 0(e). And finally, these fluctu¬ 
ations are assumed to have much larger scales along the 
magnetic field direction than across it: k\\/k± ~ 0(e), 
where fc|| and k± are the wavevectors along and across 
the magnetic held, respectively. An important conse¬ 
quence of this assumption is that the perpendicular bulk 
speed V± is mostly dominated by drifts, which have to 
be smaller than the ion thermal speed Vth,i according to 
e = V±/vth i 1. This fact will be a source of differences 
between the GK and kinetic simulations that do not have 
such a restriction. 

Simulations of magnetic reconnection with codes based 
in the GK theory have shown to drastically reduce the 
computation time, which allows the use of more realistic 
numerical parameters (e.g.: mass ratio), and at the same 
time, the possibility to capture many of the desired ki¬ 
netic effects, such as finite ion Larmor radius ones. 10-13 
However, one of the drawbacks of the GK simulations 
of magnetic reconnection is that the initial setup has to 
be force free J x B = 0. The often used Harris sheet 14 
initialization is not possible in the standard gyrokinetic 
theory, since the particles that carry the current are dif¬ 
ferent from the background, and therefore the perturbed 
current is of second order in the ordering GK parameter 
e. Nevertheless, a force free initialization is a good ap¬ 
proximation to the exact kinetic equilibrium in low beta 
plasmas with strong guide field. 

Due to the previous reasons, it is of fundamental im¬ 
portance to establish the key properties of magnetic re¬ 
connection that can be properly modeled with a gyroki¬ 
netic code, instead of running computationally expensive 
fully kinetic PIC simulations. In this context, the first 
direct comparison between PIC and GK simulations of 
magnetic reconnection in the large guide field limit was 
performed just recently in Ref. 15. They found that many 
reconnection related fluctuating quantities (such as the 
thermal/magnetic pressures) obtained for different PIC 
guide field runs have the same value after they are scaled 
to b g . And this value is equal to the result given by the 
GK simulations after a proper choice of the e parameter. 
This is valid under the assumption of a constant total 
plasma /3 (to be defined in Sec. II) among different PIC 
guide field runs. The result is in agreement with the pre¬ 
dictions of a two fluid analysis of Ref. 16. In addition, 
Ref. 15 showed that the morphology of the out-of-plane 
current density exhibits a good convergence between the 
PIC simulations with sufficiently high guide field and the 
gyrokinetic runs. They also noticed that low /3 plasmas 
require higher PIC guide fields to reach convergence to¬ 
wards the GK results, in comparison to high /3 plasmas. 

However, the previous benchmark study in Ref. 15 did 
not discuss some limitations in the GK approach, that 
might make a comparison with a fully kinetic description 
of magnetic reconnection not reliable in some parameter 


regimes. In this context, our purpose is investigating the 
physical origin of the differences between these plasma 
models in the realistic regime of finite PIC guide fields 
(since PIC cannot use arbitrarily large values of guide 
field, where results should match). Thus, we aim to es¬ 
tablish the limits of applicability of the GK compared 
to PIC simulations of magnetic reconnection. Some of 
the differences come from the special way in which the 
GK plasma model enforces the force balance in the per¬ 
pendicular direction to the magnetic field to order e. In¬ 
deed, from the perpendicular GK Ampere’s law, it can 

be proven 9 ’ 1 '^ 18 
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Here, the quantities denoted as 5 are of order e. SP± is 
the (total) perturbed perpendicular pressure tensor and 
8B\\ the perturbed parallel magnetic field. This expres¬ 
sion is satisfied exactly to order e, being formally equiv¬ 
alent to the pressure equilibrium condition as given by 
the two fluid model in the strong guide field limit (to be 

discussed in Eq. (7)). The only difference is that 8P± 
can possibly include off-diagonal elements representing 
finite Larmor radius effects in the GK model, while in 
the two-fluid model the thermal pressure is only a scalar 
quantity. Note that even though the magnetic pressure 
is much larger than the thermal pressure in this force free 
scenario, their gradients or fluctuations can still be com¬ 
parable, and that is why such pressure equilibrium still 
exists even in the strong guide field limit as assumed in 
the GK approach. On the other hand, there is no such 
restriction in the PIC approach, but it converges towards 
it as the guide field increases. 

The physical consequence of the previous fact is that 
the agreement between both approaches breaks down 
mostly when secondary magnetic islands, with strongly 
unbalanced magnetic pressures, appear in the PIC sim¬ 
ulations using a (relatively) low guide field. Although 
they also appear in the PIC high guide field regime or 
the GK simulations, these secondary magnetic islands 
have different features as a result of an initial shear flow 
present in the PIC force free initialization. These sec¬ 
ondary magnetic islands are structures that appear as 
a result of the tearing mode, with a different (smaller) 
wavelength than the one imposed by the large scale ini¬ 
tial perturbation. They start at electron length scales, 
growing up to ion scales. In this sense, our definition 
does not exactly match with some previous studies, 19-21 
which limit secondary islands to those ones with electron 
length scales. Other works have stated that secondary is¬ 
lands have opposite out-of-plane current density to those 
of the primary islands. 22 Again, the structures in our 
simulations do not match these features. 

We expect the formation of these structures in our 
setup, since it is known that guide field reconnection pro¬ 
duce a “burstier” reconnection, forming more secondary 
magnetic islands than an antiparallel configuration. 23 
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Many previous works have analyzed these structures with 
several levels of details, in particular in the context of 
extended current layers, since they can modulate tempo¬ 
rally the reconnection rates (see Ref. 24 and references 
therein). For example, the authors of Ref. 25 reported a 
core magnetic field inside of secondary magnetic islands 
by means of hybrid simulations of Harris sheets with a 
weak magnetic guide field. It was attributed to a non¬ 
linear amplification mechanism between the guide field 
and the out-of-plane magnetic field generated by Hall 
currents. Later, Zhou et al. 26 analyzed the same phe¬ 
nomenon during the coalescence of secondary magnetic 
islands, with 2D PIC simulations for a similar setup. In 
our case, we will show that a similar consequence takes 
place but as a result of a shear flow initialization in the 
PIC low guide field regime. 

Shear flows have been analyzed exhaustively for their 
importance in the development of magnetic reconnection. 
It is known that reconnection rates decrease with faster 
shear flows parallel to the reconnected magnetic field, ac¬ 
cording to the scaling law derived in Ref. 27 under a Hall- 
MHD model without guide field (basically, due to less effi¬ 
cient outflows from the X point). A kinetic approach and 
PIC simulations have confirmed this prediction, but have 
also shown that tearing mode (associated with the for¬ 
mation of magnetic islands) is never completely stabilized 
for thin CS, even with super-Alfvenic flows. 23 A two-fluid 
analytical study 29 of collisionless tearing mode driven by 
electron inertia showed that a shear flow, under certain 
conditions, can generate a symmetric out-of-plane mag¬ 
netic field. Shear flows are also closely related with the 
Kelvin-Helmholtz instability. Our PIC runs in the low 
guide field regime do not exhibit outflows with sufficient 
speed to drive the kinetic version of Kelvin-Helmholtz, 
mostly because they scale as the Alfven speed Va which 
is known to be a threshold for this instability. 30 It is also 
important to mention the work by Fermo, Drake, and 
Swisdak 31 and Liu et al. 32 , which found, by means of 2D 
PIC simulations, secondary magnetic islands generated 
by the Kelvin-Helmholtz instability in the reconnection 
outflow (similar to the finding by Loureiro, Schekochi- 
hin, and Uzdensky 33 , but in the framework of reduced 
MHD). On the other hand, the structures observed in our 
case appear only close to the main X point. However, the 
secondary magnetic islands in our PIC/GK simulations 
share some of the features seen in those previous works. 
That is why we are going to relate their appearance with 
the core magnetic field and shear flow, depending on the 
parameter regimes of our PIC/GK runs. 

The remainder of this paper is organized as follows. In 
Sec. II, we describe the simulation setup used by both 
GK and PIC codes, as well as the parameter range to 
be studied. In Sec. Ill we identified the deviations from 
the linear scaling on the guide field of some magnetic 
reconnection quantities, extending the previous compari¬ 
son work by TenBarge et al . 15 These differences between 
PIC and GK simulations, appearing for some parameter 
regimes and specific spatial and temporal locations, are 


manifested in both magnetic and thermal pressure fluctu¬ 
ations. In this paper we focus on the first ones, while the 
deviations with respect to the linear scaling in the ther¬ 
mal pressure fluctuations will be deferred to a follow-up 
paper. The remarkably high values of magnetic fluctu¬ 
ations in the PIC low guide field regime are due to the 
generation of a core magnetic field in the secondary mag¬ 
netic islands. This process is discussed in Sec. IV in the 
context of two fluid theory and violations of the pressure 
equilibrium condition. Then, in Sec. V, we analyze the 
physical mechanism of this core magnetic field as a result 
of a shear flow due to the force free initialization in the 
PIC low guide field runs. We briefly discuss the effects of 
a high /3 plasma on all these processes in Sec. VI. Finally, 
we summarize our findings in the conclusion, Sec. VII. 


II. SIMULATION SETUP 


The results of the 2.5D simulations (no variations 
allowed in the out-of-plane direction z) to be shown 
were obtained with the explicit fully kinetic PIC code 
ACRONYM 34 and the gyrokinetic GENE code. 35 For 
both PIC and gyrokinetic simulations, we used a very 
similar setup and parameter set to the ones used in 
Ref. 15. We will repeat here the more important pa¬ 
rameters to ensure a future accurate reproducibility of 
our work. The initial equilibrium is a force free dou¬ 
ble CS, with halfwidth L and magnetic field B(x) = 
B y (x)y + B z (x)z given by: 


By - Booy 


tanh 



— tanh 


^a; — 3 Lz/4^ 



B, = B a 


b 2 g + cosh 2 



+ cosh 


- 31 ^/ 4 ^ 
(3) 


with relative guide field strength b g = Bg/B^y , where 
Booy is the (asymptotic) reconnected magnetic field. This 
corresponds to a total magnetic field of constant magni¬ 
tude Bt = Baoy^Jl + bg. A change in b g is obtained by 

modifying B^y but always keeping Bt constant (note 
that Ref. 15 keeps B g constant, but Bt ~ B g if b g 1 
and so the differences are not significant in this regime). 
L x and L y are the simulation box sizes in the reconnec¬ 
tion plane x-y, respectively. Ions are initially stationary, 
while the current is carried by drifting electrons with ve¬ 
locity U e according to J = — en e U e , satisfying the Am¬ 
pere’s law V x B = /io<7. n e = rii = no are the initial 
constant electron (“e”) and proton (“i”) number densities, 
respectively. In order to accelerate the reconnection on¬ 
set, we use a perturbation in B x and B y that generates 
an X/O point at the center of the left/right CS and that 
can be derived from the vector potential 5A Z 


5A Z = SPtt~ sin 
27T 


( 2tt (y + Ly/ 4) 
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Case 

Pi 

iOpe l ^ce 

Vth,e/C 

LI pi 

Lx{= Ly) 

Low beta 

0.01 

0.8 

0.125 

2 

40-7T pi 

High beta 

1.0 

4.0 

0.25 

1 

207T pi 


Table I. Physical parameters for the two sets of runs. 


where 6P = 0.01 is the amplitude of the perturba¬ 
tion. The parameters used in this study are divided 
in two sets: “low” and “high” beta cases, being sum¬ 
marized in Table I. All the corresponding quantities 
are calculated with respect to Bt , i.e.: the ion plasma 
beta /3i = 2fj, 0 n 0 kBTi/B^, the electron/ion cyclotron 
frequency Q c {e/i} = &Bt/ m { e /i} and the electron/ion 
Larmor radius P{ e /q = ^.{e/d/^cWi}, with /lie elec¬ 
tron/ion thermal speed v th ^ e/i} = ^/2k B T {e i i} /m {e / i} . 
We also have equal temperatures for both species Xj = 
T e , implying VA = (w pe /fl ce )(u t / lie /c). All the details 
concerning normalization and the appropriate correspon¬ 
dence between GK and PIC results can be found in the 
Appendix A. 

Now, let us specify some numerical parameters. For 
all cases we use a reduced mass ratio mi/m e = 25. Both 
codes use double periodic boundary conditions (x and 
y directions). For the PIC runs, we use a grid size Ax 
such as N x = N y = 1024 cells in the low beta case (with 
Pe/Ax = 1.69), while for the high beta case is N x = N y = 
1280 cells (with p e /Ax = 4.07). The time step is chosen 
to be A tujp e = 0.03/0.12 in the low/high beta case, re¬ 
spectively, to fulfill the Courant-Friedrichs-Lewy (CFL) 
condition (for light waves) with (cAt)/Ax = 0.5 < 1. 
Thousand particles per cell are used in both cases for 
each species. A second order interpolation scheme, also 
known as TSC shape function (Triangular Shaped Cloud) 
was used to reduce numerical noise without having to in¬ 
crease significantly the macroparticle number. Finally, 
no current smoothing was used. For the GK runs, the 
spatial grid is N x = N y = 1024 for both cases, while 
the parallel/perpendicular velocity grid is chosen to be 
L v = 3vth,i , = SksTi/B t, with 32x20 points in the 
space (u||,/r). p is the (adiabatic invariant) magnetic mo¬ 
ment. In the GK simulations, the initial noise level and 
spectrum were chosen to match with the corresponding 
one in the PIC runs. Then, this noise acts as an addi¬ 
tional perturbation on top of the one described in Eq. (4) . 

Finally, some comments about the computational per¬ 
formance of the runs with our codes under the different 
parameter regimes of this study can be found in the Ap¬ 
pendix C. 


III. GLOBAL EVOLUTION AND DEVIATIONS FROM 
LINEAR SCALING 

In this section, we first show the global evolution of the 
GK/PIC simulations comparing with the results obtained 
by TenBarge et al. 15 Then, we explain and quantify the 


linear scaling on the guide field as predicted by the two 
fluid model of Ref. 16, in order to identify the deviations 
from it and the key open problems that will be addressed 
in this and an upcoming paper. 


A. Global evolution 

With our independent set of codes, the PIC code 
ACRONYM and the GK code GENE, we obtained sim¬ 
ilar values for the reconnection rates (see Fig. 1) as the 
ones reported by the previous comparison work, Ref. 15. 
First, the normalized reconnection rates have the same 
value for each high/low /3 case for different guide fields 
(although they decrease with b g when measured without 
normalization, consequence of keeping constant the total 
plasma /3). The constancy of normalized reconnection 
rates in the strong guide field limit agrees with some pre¬ 
vious studies, even when the total plasma /3 changes 32 
(and correspondingly it was expected to have different 
reconnection regimes). Note that is not valid for Harris 
CS in the very low guide field regime, where it is known 
that the normalized reconnection rates are reduced with 
increasing b g (see, e.g., Refs. 36-38), as a result of the 
enhanced incompressibility of the plasma. 

Second, we also found smaller reconnection rates in 
the high beta case compared to the low beta one, in 
agreement with previous studies of gyrokinetic simula¬ 
tions of magnetic reconnection. 39,40 However, we can no¬ 
tice in Fig. 1 that the reconnection values for both low 
and high plasma (3 cases are still a substantial fraction of 
(d’3//dt)/'i'jv ~ 0.1. Therefore, fast magnetic reconnec¬ 
tion takes place in both cases, even though the low beta 
runs are in a regime where there should be not dispersive 
waves that may explain these rates, according to the two 
fluid analysis in Ref. 41 (see also Ref. 38). This contra¬ 
diction was already addressed in Ref. 15 and also other 
recent works. 32,42,43 Although there is no a full definitive 
answer to this controversy yet, all these evidences seem 
to indicate that fast dispersive waves, such as whistler 
and kinetic Alfven waves, seem not to play the essential 
role to explain fast magnetic reconnection in collisionless 
plasmas. Since the purpose of this paper is different, we 
are not going to discuss further that issue in this work. 

It is worth to mention that some details in the evolu¬ 
tion of the system are not exactly the same as in Ref. 15. 
For instance, in the low beta case (Fig. 1(a)), although 
the peak in the reconnection rate is reached more or less 
at the same time, we observed that the CS is more prone 
to the development of secondary islands than reported 
there. They start to appear just after the reconnection 
peak ( t > 4Qta) instead of t > 75 ta- Probably, this is 
due to the different initial noise level in the PIC runs, 
since it is random and dependent on the specifics of the 
algorithms implementation. Therefore, the locations and 
timing of the secondary magnetic islands are expected to 
be different between the results given by the VPIC (as 
used in Ref. 15) and ACRONYM codes (and also with 
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Figure 1. Comparison of reconnection 
rates dip/dt for the left CS among differ¬ 
ent PIC guide field and GK runs. This 
quantity is calculated as the difference 
in the out-of-plane vector potential A z 
between the X and O points, a) Low 
beta Pi = 0.01 case, b) High beta 
Pi = 1.0 case. 


the GK code GENE). 

All the results and plots to be shown from now on 
are based in the low beta case. The slightly different 
conclusions for the high beta case will be briefly analyzed 
only in Sec. VI. 


B. Parity/symmetry of magnetic reconnection quantities and 
linear scaling 


As indicated in Ref. 15, based in the two fluid anal¬ 
ysis of Ref. 16, the thermal pressure and magnetic field 
fluctuations should display antisymmetric (odd-parity) 
structures in the separatrices in the strong guide field 
regime b g 1, assuming small enough fluctuations and 
asymptotic plasma beta p y = 2yi 0 nokBTi/B/^y > 1 (al¬ 
though a later study 29 showed that this assumption is not 
necessary for the appearance of an odd-parity magnetic 
field structure). In addition, both quantities should scale 
inversely with the guide field b g in the following way (see 
Eqs. (17) and (19) of Ref. 16): 

SPth di 1 

,,- r j-, (5) 

-*m,0 ^x Og 

SB, SP th dj pj 

B g B'^/y, o l x b g ' 


where l x is the typical length scale of variation of these 
quantities across the CS, away enough from the X points. 
The fluctuations <5 are calculated by subtracting the (ini¬ 
tial) equilibrium quantities B z (t = 0) « B g (for b g 1) 
and Pth{t = 0) := i + riiksTi x (for brevity, we 

omit the subscript _L in the thermal pressure). Note that 
in deriving the previous expressions, the pressure equi¬ 
librium condition in the limit of strong guide field has 
been used 16 (also obtained in the framework of reduced 
MHD), formally equivalent (by assuming a scalar pres¬ 
sure) to the perpendicular force balance of the more gen¬ 
eral GK equations as given in Eq. (1) (satisfied exactly 
to order e). This can be written in the following way by 
combining the previous expressions for 8P t h and 5B Z : 


SB Z _ 5P th 

B g Pl Pt h ,o 


(7) 


Note that Eqs. (5) and (6) predict that the fluctua¬ 
tions SPth and 5B Z are proportional to 1 /b g providing 


> 1 and l x /di constant for different guide fields, which 
is valid recurring to the estimate of order of magnitude 
lx/di ~ yf]3i given in Ref. 16. The last assumption re¬ 
quires additional evidence in order to be applied to our 
simulations, and that is why we proved this claim via 
the method explained in detail in the Appendix B. Us¬ 
ing the suitable normalizations explained in Appendix A, 
we show the estimates for SB Z and SPth in Fig. 2 and 
Fig. 3, respectively, for a time shortly after the recon¬ 
nection peak. These and all the other figures from the 
PIC runs shown in this paper, unless stated otherwise, 
have been averaged over t = 0.5 ta to reduce the effects 
of the numerical noise. In addition, the color scheme 
is scaled between ± the mean plus 3.5 standard devia¬ 
tions of the plotted quantity, a representative maximum 
as explained in the discussion of Fig. 5. From Figs. 2 
and 3 we confirm the result already found in Ref. 15 and 
predicted by the previously sketched two fluid model: 16 
the convergence of the PIC results towards the GK ones 
in the limit of strong guide field, in both odd symmetry 
and scaling with the guide field. As discussed by these 
authors, this is valid only in the region close to the sepa¬ 
ratrices. However, we can immediately notice a fact not 
discussed in the previous comparison work: not only the 
symmetry between both separatrices is broken in the low 
guide field regime, but also the appearance of an addi¬ 
tional magnetic field in the secondary magnetic islands 
not visible in GK or in the limit of strong PIC guide field. 
This is the main difference and new result to be the focus 
of the remainder of this work. For brevity, other differ¬ 
ences between PIC/GK will be addressed in a follow-up 
paper. 


IV. CORE MAGNETIC FIELD AND PRESSURE 
EQUILIBRIUM CONDITION 


The core magnetic field in the secondary magnetic is¬ 
lands, as well as in the y boundaries, can be understood 
from several points of view. In this section, we will de¬ 
scribe it in the framework of deviations from the two 
fluid model sketched in Sec. Ill B, based on the pressure 
equilibrium condition. 
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Figure 2. Contour plots of the scaled fluctuations in the (perpendicular) thermal pressure T5Pth/Pth,o f° r different PIC guide 
fields and GK runs, at a time t/rA = 50. T is for b g ^ef = 10 in Eq. (A3). (The same is valid for all the remaining figures of 
this paper in the low beta case.) a) PIC b g = 5, b) PIC b g = 10, c) PIC b g = 20, d) PIC b g = 30, e) PIC b g = 50, f) GK. 
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Figure 3. Same as Fig. 2 but for the scaled fluctuations in the out-of-plane magnetic field TSB z /Bt . a) PIC b g = 5, b) PIC 
b g = 10, c) PIC b g = 20, d) PIC b g = 30, e) PIC b g = 50, f) GK. 


A. Pressure equilibrium condition in the strong guide field 
limit in GK/PIC 


The core magnetic field in the secondary magnetic is¬ 
lands and y boundaries is unbalanced, in the sense that 


the thermal pressure does not decrease at these locations 
(compare Figs. 2 with 3). This is a violation of the pres¬ 
sure equilibrium condition in the strong guide field limit 
as given in Eq. (1) or Eq. (7). Since it is satisfied ex¬ 
actly to order e in GK, it is built-in in the equations that 
the GENE code solves, while PIC codes allow large de- 
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viations from it for finite guide fields, since they solve 
the full collisionless Vlasov-Maxwell system of equations. 
Note, in particular, that Eq. (7) represents a straight line 
with slope —/ 3 ,: in the plane SB Z - 5Pth- Therefore, a con¬ 
venient way to display this relation and deviations is by 
means of a 2D (frequency) histogram of these fluctuat¬ 
ing quantities, with results shown in Fig. 4(middle row). 
These plots were generated by selecting the interesting 
region close enough to the center of the CS with J z above 
10% of its initial value, indicated in Fig. 4(top row). In 
order to show the locations of deviations in the pressure 
equilibrium condition, we also plot in Fig. 4 (bottom row) 
the corresponding fluctuations in the total pressure 

^-ftotal _ Pth T Pmag -. 1 /g\ 

Pth,0 Pthfi 2 ^ ’ 

with Pmag = -B 2 /(2/xo). Fig. 4(middle row) shows that 
most of the locations in the CS are along the line —pi 
for the case PIC b g = 20, with an increasing spread in 
the lower b g regime. The corresponding GK results fol¬ 
low very accurately that line. Note a distinctive, pres¬ 
sure unbalanced, “bump” in the region SPth ~ 0 with 
5B Z > 0, being very noticeably for the case PIC b g = 5. 
Fig. 4 (bottom row) shows that these regions are mostly in 
the secondary magnetic islands and in the y boundaries. 
This excess of total pressure generates a net force towards 
the exterior of the magnetic islands, leading to an expan¬ 
sion in reconnection time scales (~ ta). <5-Ptotai increases 
going to the lower guide field regime, even though this 
quantity has been linearly scaled to b g . That is to be 
expected, since the high fluctuation level in these cases 
breaks down the small e assumption in which Eqs. (5),(6) 
or (7) (and also Eq. (1)) are based. 

Note that the histograms for low b g show a strongly 
asymmetric distribution in comparison with GK or PIC 
high b g regime (see Fig. 4(middle row)): there are more 
points located in the right bottom quadrant (SPth > 0 
and SB Z < 0) than in the left upper quadrant (SPth. < 0 
and 5B Z > 0). This is an indication of (and an efficient 
way to quantify) the asymmetry in the separatrices for 
the PIC low guide field regime (see the respective con¬ 
tour plots Figs. 2 and 3). We will explain the reasons in 
Sec. VA. 

We can summarize these findings by saying that, al¬ 
though the magnetic field fluctuations SB Z predicted by 
the GK simulations can be comparable to the PIC ones 
in the low guide field regime ( b g = 5 and 10) close to 
the separatrices, this is not true inside of the secondary 
magnetic islands or the periodic y boundaries, due to the 
deviation in the pressure equilibrium Eq. (1). 

B. Time evolution of magnetic and thermal fluctuations 

The purpose of this subsection is to analyze the dynam¬ 
ical evolution of the thermal and magnetic fluctuations, 
complementing the work by TenBarge et al. 15 . This al¬ 
low us to determine the physical mechanisms producing 


the deviations in the pressure equilibrium condition, as 
well as the asymmetric separatrices in the PIC low guide 
field regime. We quantify this by means of tracking the 
“maximum” of SPth and 8B Z in the entire region shown 
in Figs. 2 and 3. For the PIC simulations, in order to de¬ 
crease the effects of the numerical noise, this “maximum” 
is chosen as equal to the mean plus 3.5 standard devia¬ 
tions. The absolute maximum is not a good choice since 
is very prone to outlier values. In addition, the initial 
value of the respective fluctuating quantity is subtracted 
(note that it is not only the initial value Pth,o and B, 0 ), 
because: (1) a zero offset is required by GK since the 
fluctuating quantities are initially zero and (2) it mostly 
measures the numerical PIC noise, being enhanced for 
higher b g . The results are shown in Fig. 5. 

First of all, Fig. 5(a) shows that the maximum of SPth 
is reached around t ~ 50ta for all the cases, after the re¬ 
connection peak time t ~ 40ta when also the secondary 
magnetic islands start to form. 5B Z reaches maximum 
values even later. This in an indication of an additional 
mechanism generating 5B Z , different from the one pro¬ 
duced by the Hall currents due to the reconnection pro¬ 
cess itself (close to the X point, in the separatrices), and 
deeply in the non-linear phase. This is also the justifica¬ 
tion for showing the contour plots at a time t = 50ta in 
most of the figures of this paper. 

By means of the results shown in Figs. 5(a) and 5(b), 
we also confirm (complementing the work of Ref. 15) the 
convergence of the time evolution of both PIC TSPth and 
TSB Z toward the GK curves in the strong guide field 
limit b g > 20. As can be expected, the GK values follow 
a similar trend as the strongest PIC guide field ( b g = 
50). However, there are important deviations in some 
curves, noticeable earlier for smaller b g . For example, 
the PIC run b g = 5 shows deviations from the GK curve 
of TSB Z already in t > 20ta, while b g = 30 only after 
t > 45ta- Note that the deviations in TSPth are smaller 
than those of T5B Z . The (interesting) physical reasons of 
the differences for TSPth will be addressed in a follow-up 
paper. Here we explain and analyze with more detail the 
differences in TSB Z shown in Fig. 5(b). 

Although useful, Fig. 5 cannot provide us with infor¬ 
mation about the relation of these maxima with the pres¬ 
sure equilibrium condition Eq. (7). For this purpose, in 
Fig. 6 we show the 2D histograms relating these fluctu¬ 
ations in the plane SPth — SB Z for three characteristic 
times during the evolution of the case PIC b g = 5. The 
asymmetry along the pressure equilibrium line located in 
the right bottom part (SP t h > 0 with SB, < 0) starts 
well before the reconnection peak time, indicating a pro¬ 
cess that is active since the start. Because the points 
tracked in Fig. 6(a) (for t = 20 ta) are mostly located 
in the separatrices, we will see in Sec. V that it is as¬ 
sociated with the initial shear flow. On the other hand, 
the “bump” SPth ~ 0 with SB, > 0 (indicating violation 
of the pressure equilibrium condition and core magnetic 
field) starts to develop just before of the reconnection 
peak time (t ~ 30ta, Fig. 6(b)). These points are lo- 






rSPjota/Po 



rsP Tota/ p o r5P,o, a /P 0 GK 5P To|a /P 0 


Figure 4. Quantities showing the method and locations of the deviations in the pressure equilibrium condition for different PIC 
guide fields and GK runs, at a time t/rA — 50. Guide field increases to the right: ax) PIC b g = 5, bx) PIC b g = 10, cx) PIC 
b g = 20, dx) GK, where the x row is: 

x = 1 Top row: Contour plots of the out-of-plane current density J z /J z (t = 0). Magnetic held lines are shown in black contour 
lines. The region inside of the green contour satisfies J z /J z (t = 0) > 0.1. 

x = 2 Middle row : 2D (frequency) histograms with the correlation between the magnetic and thermal fluctuations TSB z /Bt 
and V8P t h/Pth,o- The diagonal black straight line is the pressure equilibrium condition Eq. (7). 

x = 3 Bottom row: Contour plots for the scaled total pressure TSP to tai/Pth,o- Note that dPtotai = 0 to machine precision in 
GK (d3). 
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t/T. 
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Figure 5. Time history of the “max¬ 
imum” value of (a): scaled thermal 
pressure YSPth/Pth,o and (b) magnetic 
fluctuations V8B z /Bt for different PIC 
guide fields and GK runs. See text for 
details about the calculation method. 


cated inside of the secondary magnetic islands or the y 
boundaries. Therefore, and different from the mecha¬ 
nism leading to the asymmetric separatrices, this is an 
indication of being caused via a process that needs to be 
build up during the evolution of reconnection. Later, in 
Fig. 6(c) for t ~ 40ta, the points spread even further 
in the “bump region”, indicating the shift of the largest 
deviations of the pressure equilibrium condition from the 
separatrices to the “bump”. In Sec. V we will see that the 
mechanism is a combined effect of both shear flow and 
magnetic islands generated via magnetic reconnection. 


GK simulations predict similar evolution for both mag¬ 
netic and thermal pressures especially during the linear 
phase of reconnection. The formation of secondary mag¬ 
netic islands breaks this similarity, producing deviations 
especially in the magnetic fluctuations 8B Z , that reach 
maxima for times much later than the reconnection peak 
time. 


We can summarize this section by saying that PIC and 
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Figure 6. 2D (frequency) histograms with the correlation between the magnetic and thermal fluctuations Y8B z /Bt and 
Y8Pth/Pth,o for the case PIC b g = 5, at the times: a) t = 20ta, b) t = 30ta, c ) t = 40ta. See Fig. 4(a2) for the calculation 
method. 


V. CORE MAGNETIC FIELD AND SHEAR FLOW 

In this section we describe the physical mechanism 
that leads to the generation of core magnetic field in the 
PIC low guide field simulations, complementing Sec. IV. 
In this point, it is important to mention some previous 
works that have found this feature, such as Ref. 26. They 
reported the generation of core magnetic field during the 
coalescence of secondary magnetic islands, as result of 
the Hall effect that twists magnetic field lines, plus flux 
transport with the associated pile up of the out-of-plane 
magnetic field. In our case, the mechanism is more re¬ 
lated to the first one but due to a different reason. 


A. Initial shear flow 


The core magnetic field seen in the PIC low guide field 
runs (see, e.g., Figs. 3(a) and 3(b)) is due to the devel¬ 
opment of strong in-plane currents growing from a shear 
flow present in the PIC initialization. Indeed, in addition 
to the component J z sustaining the CS (associated with 
B y (x )), the force free initialization with B z {x) implies, 
via Ampere’s law, the presence of an in-plane component 
J y , chosen to be carried only by the electrons, and given 

by 


Je.v — 


B 0 


tanh ^ tanh 


\ L 

u 2 (x - L x /4\ + 2 (x - 3L : 

cosh [ -— I cosh I -— 


1 ± 


Mo L 


V s + cosh" 2 ( 1 ^ /4 ) + cosh" 2 < X 3L * /4 


(9) 


This represents a counterstreaming (shear) flow of elec¬ 
trons along the center of each CS (see Fig. 7(b)), with 
maximum value V e y = —J e y /(en 0 ) oc B^y. Thus, due 
to the normalization used (see Appendix A), its magni¬ 
tude will decrease as 1 jb g (see Fig. 7(a)), while the associ¬ 
ated kinetic energy of the shear flow has an even stronger 
dependence (oc 1 /bg). Therefore, the shear flow strength 
in the PIC cases converges very quickly to the GK ini¬ 


tialization, where it is completely absent. The zero initial 
shear flow in this plasma model is because of the perpen¬ 
dicular GK Ampere’s law, (V x SB)± = V±SB Z = /j,oSJ± 
(leading to the pressure equilibrium condition Eq. (1)). 
Indeed, when J± is calculated from the force free particle 
distribution function, it turns out to be of second order 
in e. If taken into account, this would imply a SB Z of 
order e 2 , which is ruled out by construction from the GK 
equations. Therefore, any shear flow (in-plane SJ±) does 
not enter into the GK equations since they are of order 
e 2 . 

The initial shear flow is responsible for the asymmetric 
separatrices (see discussion of Figs. 4 and 6), especially 
regarding the preferential SPth > 0 over SPth < 0. This 
behavior was already pointed out by Cassak 27 , where it 
was attributed to the dynamic pressure of the shear flow, 
piling up electrons preferentially in one pair of the sepa¬ 
ratrices over the other one if strong enough, contributing 
to the increase of density, temperature and thermal pres¬ 
sure (see Fig. 2). An extended discussion about this topic 
will be given in a follow-up paper. 




Figure 7. a) Maximum initial value of the in-plane electron 
flow speed V e ,y versus b g normalized to Va and Vth,i, as well 
as the speed ratio VA/vth,i given in Eq. (A6). 
b) Initial profiles of V e , y (x)/V a across a CS centered in x = 0 
for two values of guide field, as well as the local Alfven speed 
Va(x). 


The effects of the shear flow need to be carefully con¬ 
sidered when normalized to Va (dependent on b g ) or Vth,i 
(independent on b g ). The ratio between these two speeds 
is shown in Fig. 7(a), decreasing with b g according to the 
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normalization used (see Appendix A). This speed ratio 
plays an essential role in the GK model, because of the 
perpendicular drift approximation. This is the assump¬ 
tion that the perpendicular bulk velocities V± are caused 
only by the E x B, diamagnetic VP e and other similar 
drifts. The perpendicular approximation breaks down 
when the in-plane speeds are close to the ion thermal 
speeds, such as the typically obtained in magnetic recon¬ 
nection outflows when Va ~ Vt.h,i- It can also be violated 
by the presence of waves with high enough frequency 
(possibly only captured by the fully kinetic plasma de¬ 
scription), leading to a cross-field diffusion of particles 
across the magnetic field. 44 Therefore, deviations from 
the real physical behavior of a Vlasov plasma modeled 
via PIC simulations are expected in the GK approach 
when the ratio VA/v t h,i ~ 1. This is precisely the situ¬ 
ation when comparing the GK to PIC simulations with 
b g = 5 or b g = 10 (the latter is the critical guide field for 
which VA/v t h,i ~ !)• 


B. Current/flows in secondary magnetic islands 

As a result of the initial shear flow, PIC runs with 
sufficiently low guide field build up a net vortical current 
inside of the secondary magnetic islands and the periodic 
y boundaries, as can be seen in Fig. 8. The formation of 
secondary magnetic islands wraps up magnetic field lines 
around them, deflecting the electron shear flow in the 
same direction. Therefore, a net out-of-plane magnetic 
field is generated (see Fig. 3). Note that the direction of 
the curl of J coincides with the direction of the out-of¬ 
plane magnetic field (+2 direction points into the page). 
As expected, there is a working dynamo process in these 
places, i.e.: J • E < 0, involving a transfer of energy from 
the bulk electron motion to the magnetic field (plots not 
shown here). High guide field PIC or GK runs do not 
show this effect. (Negative values of J ■ E are seen only 
near the outflows, due to the bulk motion of the plasma.) 
A comparison of this process to the dissipation J ■ E > 0 
close to the X points will be given in a follow-up paper. 

Two other consequences of the shear flow in reconnec¬ 
tion found in previous works are worth to mention in this 
point. First, the tilting of magnetic islands and genera¬ 
tion of concentric vortical flows inside of them as result 
of sub-Alfvenic flows (in agreement with our parameter 
regime) has been seen in both 2D MHD and Hall-MHD 
simulations (see Ref. 45 and references therein). Second, 
the increasing symmetry of the core magnetic field for 
low PIC b g can also be explained due to the shear flow, 
as investigated via a two fluid analysis of tearing mode 
with shear flow of Ref. 29. But there is also opposite ev¬ 
idence to this claim: Ref. 25 found, by means of hybrid 
simulations of Harris sheets, that the “S” shape of SB Z 
in the secondary magnetic islands is only consequence of 
guide field reconnection (nothing to do with shear flow). 
This might be due to the small guide field regime ana¬ 


lyzed: b g < 1, not being correct to be applied for our 
case. 

The current inside of the magnetic islands, significant 
only in the low guide field regime, is produced due to 
the Hall effect: a decoupling of motion between electrons 
and ions, as shown in Fig. 9. For the lowest PIC guide 
field run b g = 5, the ions follow the outflow due to re¬ 
connection from the X point (Fig. 9(bl) ), but the elec¬ 
trons keep their initial shear flow and are only weakly 
deflected (Fig. 9(al)) by that reconnection outflow, fol¬ 
lowing a vortical flow pattern inside of the magnetic is¬ 
lands. This characteristic (antisymmetric) flow pattern 
is barely visible for a higher guide field of b g = 20 (see 
Figs. 9(a2)-(b2)) and totally absent for the GK run (see 
Figs. 9(a3)-(b3)), where the internal flow has a symmet¬ 
ric structure following the features of the reconnection 
outflow. 

It is important to mention that because the generation 
of magnetic field is due to a Hall effect, their effects will 
be stronger when the CS is on the order of/thinner than 
Ps = \/ks{T e + Ti)/mi/Cl ci = pi = 0.1 dj, the sound 
Larmor radius (for T) = T e ), which applies very well to 
our case (L = 2 pi). Therefore, in thicker CS these ef¬ 
fects should be significantly reduced and an agreement 
between PIC and GK simulations of magnetic reconnec¬ 
tion will be more easily reached. Note that the halfwidth 
is also on the order of electron skin depth, since pi = d e , 
and so the electron inertial effects are important: 46 elec¬ 
trons are also unmagnetized (not fulfilling the frozen-in 
condition) for distances l < L/2. 

In Fig. 9, we can estimate that the speed of the elec¬ 
tron outflows from the main X point is a few times the 
asymptotic Alfven speed V etV ~ 2.2Va- It varies weakly 
with the guide field in the PIC runs. On the other hand, 
the ion outflow speeds only reach sub-Alfvenic values 
Vi jV ~ 0.8Va, and are much more constant among dif¬ 
ferent PIC guide fields and the GK runs. These values 
agree, to a first order, with a two fluid theory of mag¬ 
netic reconnection by Shay et al. 4 ‘ : the outflow speeds 
from the X point should be of the order of the in-plane 
Alfven speed V a for io ns and in-plane electron Alfven 
speed VAe = \Jmi/Tn e V a for electrons. But the initial 
shear flow is of course strongly dependent on the guide 
field (see Fig. 7). The combination of both effects deter¬ 
mines the critical PIC guide field for which the shear flow 
can generate the current that builds up the core magnetic 
field (b g < 20). 

C. Time evolution of the current and magnetic field 
generation 

Since in Sec. VB we showed that the electron/ion out¬ 
flows do not depend strongly on b g when normalized to 
Va, the in-plane current J± oc — V e> ± will display 
similar values among different PIC b g when the same nor¬ 
malization is used, as can be seen in Fig. 8. But accord¬ 
ing to the Ampere’s law, the generation of 5B Z depends 
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Figure 8. Vector plot of the in-plane current J± = J x x + J v y for two cases of PIC guide fields a) b g = 5 and b) b g = 20, at a 
time t = 50t/i. Color coded is its magnitude \J±\/J z (t = 0), with J z (t = 0) given in Eq. (A2). 
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Figure 9. Vector plot of the in-plane bulk velocity 14/,, _l for different PIC guide fields and GK runs, at a time t/7/4 = 50. Color 
coded is the magnitude |14/,,y/14i|- Guide field increases to the right: ax) PIC b g = 5, bx) PIC b g = 20, cx) GK, where x = 1 
(top row) is for the electron V4,_l and x = 1 (bottom row) is for the ion 14,x- 


on the unnormalized J±, changing with b g according to 
Eq. (A6). It can be estimated by approximating the curl 
V x B by the gradient scale length 1/A L: 

f '-f 0?) ^ 

Since B g practically does not change for different PIC 
guide fields b g (result of choosing an invariant Bt , recall 
Sec. II), by defining the constant A = / i 0 pi/B g we in¬ 
fer that (in order of magnitude), SB Z /B g ~ A J± when 
A L ~ pi , not dependent on the guide field. Note that 
this estimate relies on the fact that pi is approximately 
constant during the evolution, which in turn depends on 
the fact that T) does not have to change too much. We 


checked that the ion temperature does not increase more 
than 30% of its initial value throughout the evolution of 
the system, and therefore we can safely take pi and so A 
constant for our purposes. This is equivalent to a propor¬ 
tionality between SB Z and J j_, a relation always satisfied 
in a force free equilibrium. This is a valid first order ap¬ 
proximation since the magnetic islands grow at ion time 
scales ( 7 / 4 ), and therefore, at electron time-scales that re¬ 
lation can be satisfied quasi-adiabatically. On the other 
hand, in our runs AL can be estimated as the size across 
the x direction of either the secondary magnetic islands 
close to the X point (~ 10 pi), or the magnetic island at 
the y boundaries (~ 18/?^). Thus, the time history of the 
maximum value of both components of A Jj_ is shown in 
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Fig. 10. 

The presence of the initial J y due to the shear flow can 
be seen in Fig. 10(b), increasingly important for lower 
PIC guide fields. But the in-plane component J y , and to 
a lesser extent J x , start to grow just after the formation of 
secondary magnetic islands. There is a good agreement 
(in order of magnitude) of the maximum values that A J y 
reach with the corresponding maximum values of SB Z 
(see Fig. 5(b)) for different guide fields, justifying empir¬ 
ically our estimation Eq. (10). Note that the factor T 
should be removed in the latter for a proper comparison 
(e.g., for b g = 5, the maximum is 5B Z /B g ~ 0.034). It is 
also important to remark that the maximum values of J y 
are reached in the boundaries of the locations where SB Z 
peaks: around the secondary magnetic islands and the 
y boundaries (neglecting the region in the separatrices 
with almost zero curl). The most important conclusion 
that can be inferred from Fig. 10 is that the maximum 
values of J y , and so 5B Z , become negligible in the PIC 
higher guide field limit when measured in absolute units. 
Therefore, magnetic field generation is only effective in 
the PIC low guide field regime. 

In principle, another possibility for the core-magnetic 
field generation that it is interesting to analyze is due 
to a pinch effect and the associated magnetic flux com¬ 
pression inside of the magnetic islands (via conservation 
of magnetic flux/frozen-in condition). We checked that 
the divergence V • v± for both electrons and ions has a 
complex spatial structure and mostly odd symmetry in¬ 
side and around the islands, not correlating well with the 
properties of SB Z (plots not shown here). Moreover, we 
checked that the absolute value of V • v± (normalized to 
the natural time scale ify 1 ), is reduced for higher guide 
fields, implying a more incompressible plasma. This is 
consequence of the fact that the lowest order perpendic¬ 
ular drifts are incompressible, a condition satisfied better 
in this regime. However, V • v± decays much faster with 
the guide field as the core-field: already for b g = 20, this 
quantity is comparable with the noise level but, on the 
other hand, there is still significant SB Z in the islands for 
this guide field case. Therefore, a pinch effect V • vj_ < 0 
is not likely to be responsible for the generation of SB Z 
in the magnetic islands. 

All these are indications that the generation of the core 
magnetic field SB Z is due to a combination of the effects of 
the shear flow and the formation of magnetic islands, be¬ 
coming increasingly important for lower PIC guide fields, 
and that is why it should not be taken into account when 
comparing PIC with GK simulations of magnetic recon¬ 
nection. 


D. Effects of counterstreaming outflows in core magnetic 
field generation 

We already mentioned that in the PIC low guide 
regime, besides of the secondary magnetic islands, a core 
magnetic field is also generated in the y boundaries. This 


is due to a combination of two mechanisms. The main 
one is because of the compression by the colliding out¬ 
flows generated from reconnection in the main X point 
and the periodic boundary conditions (equivalent to a 
configuration with multiple X points). This effect is 
stronger for PIC low guide field regime since the electron 
outflows are faster (in absolute units, since they have 
same values in units of Va). Note that this process could 
be avoided by choosing longer boxes. 25 The second mech¬ 
anism contributing to the core magnetic field generation 
is the same as for the secondary magnetic islands, since 
the O point of the primary magnetic island is located 
precisely at the y boundaries. 

E. Influence of shear flow in reconnection 

As we can see in Fig. 1, the reconnection rates are 
smaller for the lower guide field cases b g = 5 and 10 com¬ 
pared to the PIC high guide field regime or the GK runs. 
This can be understood in terms of three different rea¬ 
sons. First, the outflows from the X points should be less 
efficient (slower) because of the additional magnetic pres¬ 
sure due to SB Z in the secondary magnetic islands (and y 
boundaries), after being normalized to the Alfven speed. 
This process tends to inhibit the CS thinning. Note 
that this additional out-of-plane guide field has higher 
relative importance with respect to the asymptotic mag¬ 
netic field precisely in cases with low PIC b g (see, e.g., 
Fig. 5). Second, as indicated in Ref. 27 with the Hall- 
MHD model without guide field, the magnetic tension in 
the reconnected magnetic field lines should be released by 
the shear flow, reducing the outflow speed produced by 
them and thus decreasing reconnection rates by a factor 
of (l — Ve.y/V‘ 4 )■ Similarly, numerical solutions of a ki¬ 
netic dispersion relation and 2D PIC simulations for thin 
CS 28 demonstrated a reduction in the growth rate of the 
collisionless tearing mode for increasing shear flows, the 
instability associated with the onset of magnetic recon¬ 
nection. And third, part of the available magnetic energy 
for reconnection that should be converted into particle 
energy is given back when the Hall currents are formed 
in the secondary magnetic islands. Overall, these reasons 
suggest that reconnection rates might be reduced in the 
PIC low guide field regime, if the total plasma /? is kept 
constant. 


VI. FINITE PLASMA BETA EFFECTS 

Let us finally discuss the high beta case with fa = 1.0. 
Although the basic phenomenology of magnetic reconnec¬ 
tion and secondary magnetic islands is similar, in general, 
the agreement is better between different guide field PIC 
runs in the range = 1 —>■ 10 and the corresponding 
GK results. This can be seen in the reconnection rates 
of Fig. 1(b), as well as in the time evolution of the scaled 
magnetic fluctuations T5B Z shown in Fig. 11(b). The 
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Figure 10. Time history of the maxi¬ 
mum value of the in-plane current den¬ 
sity for different PIC guide fields cases, 
a) A J x . b) A J y . Note the absolute units 
(no normalization depending on b g like 
in Fig. 8). These values have to be com¬ 
pared with the ones for 8 B Z /B g shown 
in Fig. 5. 


time evolution of the thermal pressure fluctuations T8Pth 
in Fig. 11(a) displays a different behavior between the re¬ 
sults given by both codes. As we already pointed out, the 
physical origin of these differences will be analyzed in a 
follow-up paper. 

There are at least two related reasons for the better 
agreement. The first one is because the fluctuations 
SPth/Pth, o are predicted to be smaller by one order of 
magnitude compared to the low beta case, since they are 
proportional to oc 1 /y/Wi ( see Eqs. (5) and (6)). In addi¬ 
tion, the GK ordering parameter e will also be reduced 
(see discussion of Eq. (A5)), from e = 4 for b g = 5 in 
the low beta case to e = 0.2 for b g = 1 in the high beta 
case. This improves the validity of the predictions of the 
GK compared to the PIC model in this parameter regime, 
since the largest deviations from the pressure equilibrium 
condition Eq. (7) will be reduced by the same amount. 
Also note that the maximum net force due to the total 
pressure imbalance will be much smaller than in the low 
beta case. 

Second, the range of variation of the speed ratio given 
by Eq. (A6) decreases by a small amount in comparison 
to the low beta case: VA/v t h,i = 0.1 —>• 0.7 going from 
b g = 10 —> 1. This implies that the maximum values 
of the initial shear flow (for b g = 1) are much smaller: 
max(Vyo) ~ 0.3 Va ~ 0.2 v t h,i, and so are their effects on 
the system. Moreover, reconnection rates are reduced by 
a factor of two in this regime (see Fig. 1(b)). Then, the 
maximum electron/ion outflows speeds in units of Va, 
proportional to the reconnection rates, will also be fur¬ 
ther reduced in comparison to the low beta case, becom¬ 
ing negligible when measured in units of Vth,i■ More pre¬ 
cisely, we measured maximum electron outflows speeds 
on the order of 0.75Va (see Fig.l2(3rd. row)) and ion 
outflows on the order of 0.3 Va (see Fig.l2(4th. row)), 
roughly smaller by a factor of 3 compared with the cor¬ 
responding values in the low beta case. 

This reduction in the maximum in-plane flow speeds 
has two physical consequences. For the GK model, the 
maximum deviations from the drift approximation will 
be smaller than in the low beta case (see discussion in 
Sec. VA). For the fully kinetic approach, the magnetic 
field generation due to the colliding outflows at the y 
boundaries is reduced. On the other hand, we also ob¬ 
served smaller magnetic islands in this case, implying a 
smaller core-magnetic field. This is because the genera¬ 


tion of SB Z , according to the estimate Eq. (10), is pro¬ 
portional to the length scale A L of these islands, which is 
roughly five times smaller compared to the low beta case. 
This is valid even though the Hall term and the corre¬ 
sponding decoupling of electrons and ions is facilitated 
in high plasma beta environments, implying a greater in¬ 
plane J± (twice as high compared to the low beta case). 
The net result of all these effects can be seen in the time 
evolution of T8B Z in Fig. 11(b), where the deviations of 
PIC results compared to the GK ones are significant only 
for b g = 1 and much smaller in absolute terms compared 
with the low beta case. Convergence with the GK re¬ 
sults is already reached with values b g > 3 (see, e.g., the 
second column of Fig. 12 for b g = 5). This reduction 
of core magnetic field strength in high /3 plasmas is in 
agreement with previous hybrid simulations 25 (although 
for very low guide fields b g < 1). 

Finally, it is important to mention that there is, in 
general, a higher level of numerical noise, and associated 
numerical heating, in this high beta case compared to the 
low beta one for the PIC runs. It becomes increasingly 
important for higher guide fields, reducing even faster the 
signal-to-noise ratio. This is, in part, responsible for the 
monotonically increasing thermal pressure fluctuations in 
Fig. 11(a) for later times, especially in the case b g = 10. 
This observation was already pointed out in Ref. 15, be¬ 
ing a well known consequence of the enhanced numerical 
collisions in weakly magnetized environments simulated 
by PIC codes, or equivalently, in high beta plasmas. 48 


VII. CONCLUSIONS 

We have carried out a comparison of magnetic recon¬ 
nection in the limit of strong guide field between two 
plasma models by using fully kinetic PIC with gyroki- 
netic simulations of force free current sheets. Our study 
extends the previous work of Ref. 15. We established the 
limits of applicability of the gyrokinetic approach com¬ 
pared to the fully kinetic model in the realistic regime of 
finite PIC guide fields, and the physical reasons behind 
these differences. Note that the following conclusions are 
based on sets of runs with total ion plasma /3j = 0.01 
constant for different PIC guide fields. 

First, by using an independent set of PIC and gyroki¬ 
netic codes, we found the limitations in the linear scaling 
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Figure 11. Same time histories as Fig. 5 
but for the high beta case, (a): scaled 
thermal pressure F<5 Pth/Pth,o and (b) 
magnetic fluctuations T8 B z /Bt . T is 
for fe 9 ,re/ = 5 in Eq. (A3). 


reported in Ref. 15 in both thermal and magnetic fluctua¬ 
tions. PIC simulations in the low guide field regime show 
deviations in the inverse scaling with the guide field, not 
converging properly to the gyrokinetic results. These de¬ 
viations start to be especially significant when secondary 
magnetic islands start to form, after the reconnection 
peak time. In this paper, we focus on the additional 
magnetic fluctuations revealed by the PIC low guide field 
runs, which are mostly due to macroscopic bulk plasma 
motions. The analysis of the differences in the thermal 
fluctuations—that have to do with dissipative processes, 
heating mechanisms and non thermal effects—will be ad¬ 
dressed in a follow-up paper. 

In particular, we found that the PIC low guide field 
regime shows an excess of magnetic pressure inside of the 
secondary magnetic islands (core magnetic field), reach¬ 
ing maximum values for times much later than the re¬ 
connection peak time. Correspondingly, the agreement 
between the two plasma descriptions is better before the 
formation of these structures, during the linear phase of 
magnetic reconnection. This excess of magnetic pressure 
is not compensated by a corresponding decrease in the 
thermal pressure. The reason is that gyrokinetic codes 
keep the pressure equilibrium condition to machine pre¬ 
cision because perpendicular force balance is enforced in 
the GK equations, while PIC codes allow large devia¬ 
tions from it since they solve the full collisionless Vlasov- 
Maxwell system of equations. Therefore, although the 
gyrokinetic results can be comparable to the PIC ones 
for relatively low guide fields ( b g = 5 and 10) in the 
separatrices close to the X points, convergence in the 
secondary magnetic islands requires much higher guide 
fields (b g > 30). 

We found that the physical mechanism that generates 
the core magnetic field (associated with a dynamo ef¬ 
fect J ■ E < 0) is due to an initial shear flow present in 
the force free PIC initialization. This shear flow, that 
also produces asymmetric separatrices and reduces re¬ 
connection rates, is negligible in the limit of strong guide 
field and absent in the corresponding gyrokinetic initial¬ 
ization. Through the Hall effect that decouples electron 
from ion motion, vortical electron flow patterns are gen¬ 
erated in the secondary magnetic islands as magnetic re¬ 
connection wraps up magnetic field lines around these 
structures, carrying the electron shear flow with them. 
This magnetic field weakens for higher PIC guide field 


runs, converging to the gyrokinetic result, since the shear 
flow is not strong enough to drive the currents capable of 
generating it (when measured in absolute units). There 
is also an out-of-plane magnetic field generation at the 
periodic y boundaries, where it is located at the O point 
of the primary magnetic island. This is not only due to 
the same previous mechanism, but also due to the com¬ 
pression by colliding electron outflows, faster in the PIC 
runs with low guide field (in absolute units). 

We also showed that the relative ratio of the electron 
outflow speed to the ion thermal speed, proportional to 
VA/vt.h,ii is higher for the PIC low guide field regime, 
reaching values close to 1 for b g ft 10. This breaks the 
perpendicular drift approximation, a critical assumption 
in which the gyrokinetic codes (and the gyrokinetic the¬ 
ory in general) are based, and thus an additional source 
of differences is expected compared to the PIC simulation 
results. 

Finally, we also analyzed the effects of a high plasma 
beta ft = 1.0 compared to our standard case ft = 0.01. 
Although the basic phenomenology is similar, we could 
notice a better agreement between the corresponding PIC 
and gyrokinetic results, reaching a relatively good con¬ 
vergence for guide fields as low as b g ~ 3. From this, 
we can conclude that an accurate comparison between 
PIC and gyrokinetic force free simulations of magnetic 
reconnection requires high plasma /3 ~ 1 (to reduce the 
fluctuation level proportional to 1/ftft), although PIC 
codes are affected by enhanced numerical heating in this 
regime. Moreover, it is necessary to have parameters 
with low ratios VA/vth,i 1 (to avoid the effects of the 
initial shear flow), a small ordering parameter e C 1 
in the gyrokinetic initialization, and reconnection rates 
should not be too high (dT/ dt)/ 4 /n ft 0.1. In the PIC 
model, this is to avoid super-Alfvenic outflows generat¬ 
ing stronger magnetic fields at the periodic boundaries, 
while in gyrokinetic model these flows might also break 
the perpendicular drift approximation. The balance of 
these related numerical and physical parameters allow to 
determine a convenient parameter regime for an accu¬ 
rate comparison of magnetic reconnection between these 
plasma models. 

In spite of all those differences, the gyrokinetic simu¬ 
lations show a development of magnetic reconnection re¬ 
markably similar to the PIC high guide field runs. This is 
of central importance due to the computational savings 
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Figure 12. Contour plots showing some quantities in the high beta case Pi = 1.0 for two PIC guide helds and GK runs, at a 
time t = 50ta- Guide field increases to the right: ax) PIC b g = 1, bx) PIC b g = 5, cx) GK. Row x = 1: scaled total pressure 
rdPtotai/ Pth., 0 - Row x = 2: scaled magnetic fluctuations F5B Z /Bt- Row x = 3: Vector plot of the in-plane electron bulk 
velocity with color coded the component V RiV /Va- Row x = 4: Vector plot of the in-plane ion bulk velocity with color coded 
the component V.m/Va- 


of the first ones. Indeed, for the low beta case, we mea¬ 
sured speed-ups by a factor of 10 3 between the gyroki- 
netic runs compared to a corresponding PIC guide field 
b g = 50 simulation. However, the computational cost of 
a PIC simulation decreases linearly with the guide field, 
in such a way that in the low guide field regime ( b g ~ 5), 
the speed-up and advantages of using gyrokinetic instead 
of PIC simulations are not so notorious. 
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Appendix A: Normalizations 

Here, we explained all the details regarding normaliza¬ 
tion and the right choice of a correspondence between the 
GK and PIC results. 

The lengths are normalized to pi and the times to the 
Alfven time ta = L/Va, where Va is the Alfven speed 
defined with respect to the reconnected magnetic field 

Va _ 1 Bpoy _ _1_ 

c c VMonom* ^ + b 2^ pe /n ce )^nJwT e ' 

(Al) 

where ui pe = ^/n e e 2 /(£om e ) is the electron plasma fre¬ 
quency. We use the previous definitions for the normal¬ 
ization of the reconnection rate: i/n = B^yV a and cur¬ 
rent density: Jn = enoYA\fWi = enoVth,i/\J 1 + b 2 . Note 
that due to the Ampere’s law, the normalized initial out- 
of-plane current density that supports the asymptotic 
magnetic field B y is given by: 

_ J z (t = 0) _ en 0 U e _ U e j 2 _ 1 

zN '~ Jn ~ Jn ~ v th>i V + (L/d^VA’ 

(A2) 


which is independent on the guide field strength. 

As we mentioned in Sec. IIIB, many of the fluctuating 
quantities scale with the guide field due to the fact that 
/3i, and so the transverse distance l x /di , are constant for 
different b g (with the last claim proved in Appendix B). 
Therefore, the PIC fluctuating quantities will have the 
same value under the previous assumption if they are 
multiplied by the factor proportional to b g 


r = 


\J l + b l 


O g ,ref 


(A3) 


where b gre f is a reference guide field, chosen to be 
b g = 10/5 in the low/high beta case, respectively. The 
value b g = 10 in the low beta case was chosen because 
the respective runs are not so dominated by numerical 
noise as for larger guide fields ( b g > 30) and so it is 
easier to compare with the noiseless GK runs. In addi¬ 
tion, b g , ref = 10 is not in the lowest end of guide fields 
like b g = 5 where other effects not captured by GK (ex¬ 
plained in the results, Secs. Ill, IV and V), dominate the 
physics of the system. In this way, using Eqs. (A3) and 


(5), the quantity that should be equal among different 
guide field PIC runs is 

J 1 + b l a / 1 + b l 

r SP th = \ - 5P th = V in SP th , (A4) 

Dg,ref J-U 


and analogously for SB Z . The factor y 1 + b 2 instead of 

b g is in order to keep the total magnetic field Bt constant 
(and not only B g ). 

On the other hand, the definition of T in Eq. (A3) also 
fixes the ordering parameter e of the GK runs, defined 
by: 


e 


1 


®oo y,normbg,ref 


(A5) 


where 6 003/inorrn = B ooy /B ooy re f is the normalized 
asymptotic magnetic field with respect to a reference 
value Root/, re/ expressed in code units. The initializa¬ 
tion in the GK runs gives 6ooy,norm = 0.05/2.5 for the 
low/high beta cases, respectively. Therefore, using the 
aforementioned PIC b g ^ re f, we have e = 2/0.04 for the 
low/high beta case, respectively. Note that even though 
in the low beta plasma regime with e > 1 the GK ordering 
formally does not hold, and in agreement with Ref. 15, 
this plasma model can still make accurate predictions, 
but only under some circumstances (clarified in the re¬ 
sults, Secs. Ill, IV and V). Moreover, the unusually high 
value of e is required to have a perturbed current strong 
enough to generate the relative large perturbed asymp¬ 
totic magnetic field for the case b g = 10 (being reduced 
for higher guide fields). 

The choice of a constant for different guide fields has 
another important consequence for the ratio of Alfven to 
ion thermal speed. Indeed, due to that assumption, the 
PIC runs will invariably have to change this parameter 
for different guide fields 


Va = 1 1 

yjl + b 2 g VA 


(A6) 


Note that v t h,i is constant for different guide fields, while 
Va scales inversely with it. This means that the ra¬ 
tio Va/V th,i increases for lower PIC guide fields: from 
VA/v t h,i = 0.2 —/ 1.96 going from b g = 50 -* 5 in the low 
beta case. 

It is also important to mention another side effect of 
the normalization used (already noticed in Ref. 15). Since 
Booy decreases for increasing b g , it becomes more com¬ 
parable to the noise level for very high b g . Then, for 
large b g the signal-to-noise ratio of all the fields can be 
very low, as can be seen in the extreme case b g = 50 in 
Fig. 2(e) and Fig. 3(e) (even with the extended tempo¬ 
ral average). Therefore, although in principle there is no 
limitation due to numerical constraints on the electron 
gyromotion for even higher guide fields (in our setup, p e 
and fl ce are constant for different guide fields), the PIC 
results in this regime are not reliable due to this unfor¬ 
tunate fact. 
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Appendix B: Transverse distance 


Here, we proved the statement written in Sec. IIIB: 
the transverse distance of the thermal and magnetic fluc¬ 
tuations scale as l x /di ~ y/Jk- This implies: (1) it is 
constant for different guide fields (requirement for the 
inverse scaling with b g ) and (2) it is higher for the high 
than the low beta set of parameters. We estimated the 
value of lx/di shown in Fig. 13 by using the plots of 6 Pth 
in Fig. 2. First, we detected the main X point by locating 
the minimum of the vector potential A z along the center 
of the CS. Then, we chose several distances to the left 
of that point in the y direction along the center of the 
CS ( l y , indicated in the x axis of Fig. 13). From each 
point we measured the transverse distance across the x 
direction, l x , from the center until the point in which 
the current density J z drops to 20% and 10% of its ini¬ 
tial peak value, as a simple way to detect the boundaries 
of the CS in these regions (averaging in both positive 
and negative x directions). These values correlate well 
enough with the approximate boundaries of the regions 
with significant values of SPth (as well as SB Z ) above 
numerical noise, besides of being independent of scaling 
arguments. The small error bars are the differences be¬ 
tween the 10% and 20% of the initial value of J z . First 



y distance from X point [pj 



Figure 13. Estimation of the transverse distance l x /di for 
both low (left panel) and high (right panel) beta case and 
different guide fields. These values were taken from plots of 
the thermal pressure SP t .h after the reconnection peak time 
t = 50 ta shown in Fig. 2 (for the low beta case). The specific 
method is explained in the text. 

of all, in Fig. 13(left) we can see that for the low beta 
case, the transverse distance is more or less constant for 
different guide fields at least for (longitudinal) distances 
l y < 15 pi away from the main X point, confirming the 
assumption stated before. As can be expected, the trans¬ 
verse distance l x increases away from the X point, conse¬ 
quence of the more open separatrices but also because the 
distinction between them and the boundary of the main 
magnetic island is more diffuse. In this region, there 
are more deviations from the previous constant trend for 
l y < 15 pi among different guide fields, but in any case, 
they are not significant, and a value l x ~ 0.3 di, cor¬ 
responding to 3 times the order of magnitude estimate 
lx/di ~ y/pi = 0.1, can be considered representative not 
too far away from the X point. Note that the deviations 


from more or less constant values are specially significant 
for the cases of higher guide field b g = 30 or b g =50, since 
the structures of the separatrices away from the X point 
and close to the O points are different from the low guide 
field cases (see Fig. 2). On the other hand, the estima¬ 
tions for the transverse distance for the high beta case in 
Fig. 13(right) show larger variations among guide fields, 
since in many cases the CS develop small secondary is¬ 
lands to the left and very close of the main X point and 
so the measurement is biased. In addition, for this high 
beta case, the enhanced level of numerical noise makes 
the detection of the transverse distance l x more inaccu¬ 
rate. For this case, it is only possible to conclude that the 
transverse distance is in between a certain range with a 
spread of A l x /di ~ 0.5. Nevertheless, close enough to the 
X point (l y < 5 di), l x = 1.3—1.5 di can still be considered 
reliable enough for a comparison (excluding the extreme 
case b g = 10). This value turns out to be practically, in 
order of magnitude, the estimate l x /di ~ y//5i =1.0. 

Therefore, from the results shown in Fig. 13, we can 
confirm that l x /di is (1) more or less constant for different 
guide fields, especially in the low beta case, and (2) its 
order of magnitude is between 1-3 times ~ yfj5i , which is 
good enough to ensure the validity of the inverse scaling 
with b g of 6Pth and SB Z . Note that if we had chosen 
to keep Bscy constant and increase B g to have a higher 
guide field effect, the total f3i would not be constant and 
so the inverse scaling, implying that a direct comparison 
between the different PIC guide fields and GK runs would 
not be possible. 


Appendix C: Computational performance of the GK/PIC 
codes 


It is worthwhile to mention the speed-up of the GK 
simulations compared to the PIC runs, the main motiva¬ 
tion of using the first numerical method over the second 
one for magnetic reconnection studies. Because of the 
choice of keeping the total plasma /3 constant for dif¬ 
ferent PIC guide fields to properly compare with the GK 
results (see Appendix A), the Alfven time in units of 
is (practically) linearly proportional to b g according to 


^~A^pe 


yjTi/T e mi/m e 

Vth,e/C 


fJl v 


(Cl) 


In the PIC runs, because the time step has to be propor¬ 
tional to for stability reasons, the latter expression 
will also be proportional to the number of time steps 
used to reach a given time measured in ta, and thus 
to the computational effort. Then, PIC high guide field 
runs are computationally more expensive than the runs 
in the low guide field regime (as well as the high beta 
cases compared to the low beta ones). In fact, for the 
low beta case, the ACRONYM PIC code used 3.38 • 10 4 
CPU core-hours to run the cases b g = 5 up to ta = 70 
(the last time shown in Fig. 1), while 3.33-10 5 CPU core¬ 
hours for the case b g = 50. A significant fraction of this 
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computational effort is spent in running time diagnostics 
for higher order momenta of the distribution function, 
which are used in the results to be shown in a follow-up 
paper. On the other hand, the single GENE GK code 
simulation with which those runs were compared used 
only 350 CPU core-hours, representing a speed-up by a 
factor of 10 3 when comparing with the largest PIC guide 
field considered b g =50. These huge computational sav¬ 
ings are an additional justification for the importance of 
a proper comparison study of GK with PIC simulations 
of guide field reconnection, one of the purposes of the 
present paper. 
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